Assessing diversity of King Crab Lithodes spp. in the south‐eastern pacific using phylogeny and molecular species delimitation methods

Abstract The purpose of this study was to test the hypothesis that the genetic diversity of commercially significant species of King Crabs (Lithodes spp.) along the south‐eastern Pacific (SEP) comprises different independent evolutionary units (IEUs) with spatially isolated distribution. Nine localities from inner and open waters along the SEP Chilean coast (39°S‐55°S) were sampled. We analyzed sequences from 173 individuals for the mitochondrial gene Cytochrome oxidase I (COX‐I), 151 individuals for the Internal Transcribed Spacer 1 (ITS) and 135 for the structural ribosomal RNA (28S). Genetic delimitation was performed through three analytical methods: ABGD, GMYC, and its Bayesian implementation, bGMYC. Bayesian phylogenetic analyses and haplotype networks were also performed. Divergence time between clades was assessed for the COX‐I marker and estimated from known evolutionary rates for this marker in other crustacean species and fossil calibration from other Anomuran species. Delimitation analyses, phylogenetic analyses, and mitochondrial haplotype networks suggested the presence of two deeply divergent mitochondrial lineages of Lithodes in the SEP, referred to as Clade1 and Clade 2. Nuclear markers showed low phylogenetic resolution and therefore were unsuitable for molecular species delimitation. Divergence time analysis of the mitochondrial lineages suggests a separation between Clades of approximately 2.3 Mya. The divergence time obtained suggested that Pliocene glaciations and deglaciations cycles could be involved in hybridization events between Lithodes IEUs at southern tip of South American coasts. The different frequencies of Lithodes haplotypes in inner and open water environments along SEP coasts could be explained by events such as the last glacial maximum or by differences in the adaptation of each clade to different environments. These findings support the necessity of evaluating the taxonomic status of Lithodes individuals found along SEP coasts under an integrative taxonomy approach or through markers with other evolution rates than those already used.

the presence of two deeply divergent mitochondrial lineages of Lithodes in the SEP, referred to as Clade1 and Clade 2. Nuclear markers showed low phylogenetic resolution and therefore were unsuitable for molecular species delimitation. Divergence time analysis of the mitochondrial lineages suggests a separation between Clades of approximately 2.3 Mya. The divergence time obtained suggested that Pliocene glaciations and deglaciations cycles could be involved in hybridization events between Lithodes IEUs at southern tip of South American coasts. The different frequencies of Lithodes haplotypes in inner and open water environments along SEP coasts could be explained by events such as the last glacial maximum or by differences in the adaptation of each clade to different environments. These findings support the necessity of evaluating the taxonomic status of Lithodes individuals found along SEP coasts under an integrative taxonomy approach or through markers with other evolution rates than those already used.

| INTRODUC TI ON
The importance of biodiversity and distribution of species surveys lies in the fact that these factors are the basic units for sustainable management and biodiversity conservation (King, 2007). This involves the identification and delimitation of species; from there, patterns of genetic subdivision of the species in the geographical space, which are crucial to understanding events at the population level, can be studied (Carvalho & Hauser, 1995;Carvalho & Nigmatullin, 1998). In the case of exploited resources, this information is of great interest to fishery managers since populations are the biological units on which management and conservation strategies are applied (Carvalho & Hauser, 1995;Waples et al., 2008).
Therefore, knowledge of the diversity and distribution of species targeted by the fishery can improve fishery management plans by better defining their operational units.
However, species delimitation can be a delicate task as not all characteristics will vary at the same time and not all criteria and methods used to assess these characteristics will result in the same species delimitation (De Queiroz, 2007). Due to the large number and incompatibility of some species concepts (De Queiroz, 2007;Mayden, 1997), a good option would be to refer to independent evolutionary units (IEUs, Jones, 2017;Lim et al., 2012) instead of species, especially if only one or a few concepts are considered to define them.
The delimitation of IEUs, which lies at the core of taxonomy and ecology, is related to the boundary between micro-and macro-evolution and determines number and boundaries of IUEs (De Queiroz, 2007). By analyzing DNA sequences, IEUs can be identified as "separately evolving metapopulation lineages" (sensu De Queiroz, 2007). Recent developments in analytical methods in phylogenetics, namely, relationships among lineages and the membership of individuals within these groups, can now be evaluated including the generalized mixed Yule-coalescent (GMYC) model (Pons et al., 2006), its Bayesian implementation (bGMYC) (Reid & Carstens, 2012), and the Automatic Barcoding Gap Discovery (ABGD) model (Puillandre et al., 2012), among others. Furthermore, the delimitation of IEUs based on only one type of genetic marker and one type of approach could be inconclusive, and the use of a multilocus analysis including different statistical approaches is thus highly recommended so as to obtain more precise and accurate results (Carstens et al., 2013;Dudgeon et al., 2012;Moore, 1995;Reid & Carstens, 2012). Southern King Crab (SKC) fishery is recognized as an activity of high social and economic importance in Chile (Bozzeda et al., 2019;Molinet et al., 2020;Nahuelhual et al., 2018).
Between 2016 and 2019, this fishery reached landing levels of 5179 mean tons annually (Servicio Nacional de Pesca y Acuicultura, SERNAPESCA, 2020); the haul is primarily destined for international markets, primarily China (Enexpro project, 2017). Until 2018, this fishery registered 7115 artisanal fishers (Subsecretaría de Pesca y Acuicultura, SUBPESCA, 2018) while many other workers are employed in the related industrial and commercial areas. The SKC fishery operates along the SEP coast from Valdivia (39°S) to Cape Horn (56°S), including channels, fjords, gulfs (inner waters areas), and also areas offshore (open waters). The fishery range coincides with the range described for SKC individuals along the SEP coast (Bozzeda K E Y W O R D S genetic diversity, Lithodes diversity, mitochondrial lineages, molecular delimitation methods, south-eastern Pacific

T A X O N O M Y C L A S S I F I C A T I O N
Evolutionary ecology; Genetics F I G U R E 1 Dorsal (left) and frontal (right) view of a male specimen of the morphospecies Lithodes santolla Molina 1782. Photographs taken during sampling in Tenaun, Chile, 2018Chile, . et al., 2019Molinet et al., 2020;Retamal & Moyano, 2010). Despite of the importance of this fishery, studies that describe the biodiversity and distribution of Lithodes IEUs in this area are scarce (e.g., Retamal, 2012). Currently, based on classical morphological characteristics, the SKC constitutes a single species, Lithodes santolla Molina, 1782 (Crustacea: Anomura, Figure 1) distributed along the SEP coast (e.g., Lovrich, 1997;Sierpe & Sanhueza, 2003) co-habiting with Lithodes turkayi Macpherson, 1988 andLithodes confundens Macpherson, 1988 at the southern tip of South America (Boschi & Gavio, 2005;Lovrich, 2014;Lovrich & Tapella, 2014). In a recent phy- Given the importance of clarifying the biodiversity composition and the spatial distribution of Lithodes throughout the SEP for both conservation and fishery management issues, in this study we aim to delimitate IEUs of Lithodes through an exhaustive sampling collection along the SEP coast. Three molecular genetic markers, one mitochondrial and two nuclear, and four different approaches are considered for the genetic delimitation analysis.
We hypothesize that IEUs of Lithodes along the SEP coast are constituted by at least two genetically distinct lineages with different spatial distribution. Irvine, CA, USA). We considered three molecular markers with different evolutionary rate, one mitochondrial marker, the Cytochrome Oxidase I (COX-I), and two nuclear markers, the Internal Transcribed Spacer 1 (ITS-1) and the structural ribosomal RNA 28S (28 S) for the large subunit (LSU). The universal primers defined by Folmer et al. (1994) were used to amplify the COX-I. The primers defined by Chu et al. (2001) were used for the ITS-1. These genetic markers were selected as they have been previously used for delimitating Lithodes species (e.g., Hall & Thatje, 2018;Noever & Glenner, 2018) and for their availability in public database GENBANK. The primers to amplify the gene 28S were designed in the Geneious version

| Phylogenetic analyses
Phylogenetic reconstructions for each molecular marker and for the concatenated data were conducted under Bayesian inference (BI) with MrBayes 3.2.7 software available on the CIPRES online platform (Miller et al., 2010). Substitution models for each marker were obtained with MEGA v10.2.5 (Kumar et al., 2018) and selected based on the Akaike Information Criterion (AIC) and CIPRES platform options: JC (Jukes & Cantor, 1969) for 28S, JC + I for ITS1 and HKY + G for COX-I. In Bayesian analysis with concatenated dataset, we set specific substitution models for each marker. Monte Carlo Markov Chains (MCMC) were carried out with 10 million generations, with samples of phylogenetic trees taken every 1000 generations. We discarded the first 25% of trees as "burn-in" and the remaining trees were used to generate a majority-rule consensus tree. Temperature parameter (final Temp value) was maintained to default value of 0.5.
To reveal the genealogical relationship among the mitochondrial haplotypes, networks were generated using HapView (Salzburger et al., 2011).

| Molecular delimitation
IEUs delimitation analysis considered the concatenated dataset and three different approaches. The ABGD (Puillandre et al., 2012) does not require a prior phylogenetic tree and based on genetic distance values, automatically finds the location of a barcode gap between candidate species or between intraspecific and interspecific diversity. The ABGD analysis was run using the web version (http://www.abi.snv.jussI EUs.fr/publi c/abgd/abgdw eb.html) with the best substitution model obtained for the concatenated dataset in MEGA v10.2.5 (Kumar et al., 2018), Kimura's two-parameter model (K2P model), and with priors that ranged from Pmin = 0.001 to Pmax = 0.1 with 10 steps. Next, the GMYC (Pons et al., 2006) single threshold model was carried out using the Species Limits by Threshold Statistics approach in the Splits and Ape packages of R program v. 4.0.5 (www.r-proje ct.org). This method is based on the approach that delimits species by adjusting ramified models of intra and inter species to a reconstructed gene tree (Reid & Carstens, 2012). Finally, the bGMYC was also carried out in R program v 4.0.5 with the bGMYC package (Reid & Carstens, 2012). This Bayesian analysis was based on samples from the posterior distribution of gene trees, thus allowing uncertainty in both the topology and branch lengths to be reflected in posterior parameter estimates (Reid & Carstens, 2012). A range of probabilities >0.95 was considered as strong evidence that the groups compared were conspecific, while a range of probabilities <0.05 strongly suggested that the groups compared were not conspecific (Reid & Carstens, 2012).
Inter-lineage genetic distance for COX-I dataset was calculated using Kimura's two-parameter model (K2P model) in MEGA v10.2.5 (Kumar et al., 2018). a relaxed molecular clock with lognormal distribution (Drummond et al., 2006) to mtDNA sequences using a substitution model of HKY + G obtained with MEGA v10.2.5 (Kumar et al., 2018) and an uncorrelated-lognormal (ucln) model of molecular evolutionary rate. A birth-death speciation prior was used to estimate branching rates in the phylogeny. Taking into account published substitution rates for COX-I in other crabs (Ketmaier et al., 2003;Schubart et al., 1998;Sotelo et al., 2009;Xu et al., 2009), we set the clock rate to 1.0E-8. Ucld.mean prior was set to gamma distribution with an initial value of 1.0E-8, scale of 1000 and Shape of 0.001.

| RE SULTS
The number of total sequences obtained for each genetic molecular marker are detailed in Table 1. Analysis of our sequences alignments revealed 40 variable sites (38 in COX-I and two in ITS-1), of which 37 were parsimony informative (35 in COX-I and two in ITS-1). In ITS-1, the two informative sites correspond to substitutions. Alignment lengths for each genetic marker are specified in Table S1.1. The sequences were deposited in GenBank (Accession number for COX-I: ON807360-ON807532; 28-S: ON868924-ON869058; ITS-1: ON869067-ON869217).

| Phylogenetic analysis
The number of sequences for each genetic molecular marker are detailed in Table 1. Alignment lengths for each genetic marker are specified in TableS1 The phylogeny evidence co-occurrence of both clades among almost all locations. However, spatial trends in the frequencies of Clade 1 and 2 were noted ( Figure 2). Individuals collected in the southernmost localities belong almost exclusively to Clade 1 TA B L E 1 Sampling localities, coordinates and number of sequences obtained for each genetic marker and for the concatenated dataset by locality. A haplotype network based on COX-I showed a 20-step mutational separation between the mitochondrial Clades 1 and 2 with 17 haplotypes identified within each clade (Figure 3). Interestingly, the most common haplotypes within each Clade were shared by individuals from different locations. The most common haplotype in Clade 1 is shared almost entirely by individuals from inner water locations (Teknion, Calbuco, Bahía Águila, and Fiordo Yendegaia). Clade 2 has two frequent haplotypes in the center of the network containing individuals from all sampled locations.

| Molecular IEUs delimitation
The ABGD analysis ( Figure 2)  The bGMYC resulted in a significant delimitation of Clades 1 and 2 as distinct genetic clusters, with most pairwise comparisons having a non-conspecificity probability (p ≤ .05). The mean genetic distance between Lithodes mitochondrial clades recovered was 4.58% (SD = 0.92).

| DISCUSS ION
Delimitation approaches tested here based on three different genealogical analytical methods (ABGD, GMYC, and bGMYC) showed similar results and consistently identify two genetic clusters of Lithodes in SEP coasts. Inter-lineages divergence values obtained with concatenated data using the COX-I gene were similar from others found among some crustacean species (e.g., Meyer et al., 2013;Naim et al., 2020). Phylogenetic and lineages delimitation analyses were consistent and revealed a high divergence between the two clades (Clade 1 and Clade 2). It is worth mentioning that COX-I mitochondrial marker was the main determinant of the divergence observed between clades. In this sense, the present results are coherent with those from Pérez- Barros et al. (2015), with the presence of two sympatric mitochondrial lineages of Lithodes inhabiting SEP coasts.
Nuclear markers did not show strong phylogenetic resolution in the individual phylogenies for each marker nor in the analysis with the concatenated dataset. This suggests that the nuclear markers here used may not be variable enough and new molecular diagnostic approach using a wide genomic screening perspective (like SNPs) could be the next step to clarify this point (Dufresnes et al., 2021;Herrera & Shank, 2016;Ortiz et al., 2021). whereas the Atlantic coast was completely ice free (Rostami et al., 2000).
Along these lines, it would be interesting to conduct studies on habitat preference and/or physiological requirements of each from the local effects of selection (Balding & Nichols, 1995;Vitalis et al., 2001).
A serious limitation to the utility of molecular markers as a practical resource for species diagnosis is the human error and uncertainty in creating and curating reference libraries (Collins & Cruickshank, 2013). One of the complexities in our analysis is that the sequences for L. santolla and L. confundens available in the GenBank database were formerly named based on the morphology and sampling location, which could lead to misunderstandings when trying to clarify the taxonomic status of cryptic lineages of The delimitation of Lithodes IEUs on SEP should be extended to a more integrative taxonomy study that could include spermatic analysis, crossbreeding experiments, cytogenetic analysis, and comparisons of internal morphological traits less exposed to environmental factors (i.e., spermio-taxonomy, gastric mills); this could help to identify the spatio-temporal limits of the genetic lineages identified in our study.
In conclusion, we report evidence of two deeply divergent mitochondrial lineages of Lithodes from Valdivia (39°S) to Cape Horn (56°S) along the SEP with different zonal and meridional frequencies. Considerations for fishery management should recognize the mitochondrial genetic diversity observed here, with a special focus on those localities where this diversity is greatest. In this sense, it would be interesting to carry out phylogeographic analyses at a finer scale with other more variable molecular markers to define population boundaries for harvest management and recognize biodiversity hotspots in order to conserve them (e.g., Red King Crab in Alaska, Grant et al., 2014). However, before any consideration, the taxonomic status of Lithodes on SEP coasts should be re-examined. The spatial complexity reported here, (i.e., different mitochondrial frequencies found in inner and open water areas with evidence of exchange of individuals among locations and latitudinal differentiation in haplotypes frequencies of each clade) underscores the urgency of understanding the evolutionary history of the Lithodes spp. across the SEP coasts.

ACK N OWLED G M ENTS
We are grateful for the samples provided by the FONDECYT project 1170507 "Dinámicas espaciales y batimétricas de Lithodes santolla (Decapoda, Lithodidae) (Molina, 1782) en canales del Sur de Chile: Bases para el manejo pesquero" and by the National Commission of Scientific and Technological Investigation of Chile through the FONDAP program research center (IDEAL): Dynamics of High-Latitude Marine Ecosystems (grant n. 15150003). We are also grateful for the contributions to scientific meetings and courses provided by the IDEAL project. We also thank the FONDECYT project 1220179.

CO N FLI C T O F I NTE R E S T
None declared.

DATA AVA I L A B I L I T Y S TAT E M E N T
The sequences of dataset were deposited in Genbank; COX-I: ON807360-ON807532; 28-S: ON868924-ON869058; ITS-1: ON869067-ON869217.